home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / PROGRAMM / CC_C / 0924.ZIP / SMFPS.LIB < prev    next >
Text File  |  1988-01-10  |  9KB  |  236 lines
  1. /* SMC float function library */
  2. #list-
  3.  
  4. /* LOGARTIHMIC AND TRIGONOMETRICAL FUNCTIONS */
  5. /* Trig    functions are all in degrees measure */
  6.  
  7. float exp (x)  float x;
  8.     { float     fint(), fabs(), chebgen();
  9.       static float e, m, n,    expcc[]    =
  10.       {  1.4569998748 E+0,     2.4876243388 E-1,   2.1446555995 E-2,
  11.          1.2357140820 E-3,     5.3453058176 E-5,   1.8506907121 E-6,
  12.          5.3411876864 E-8,     1.3215163808 E-9,   2.8613237815 E-11 };
  13.       static int v,    *ex = &e, *mx =    &m,
  14.              ns    = sizeof (expcc) / sizeof (float);
  15.       if (fabs (x) >= 400.)    return (x < 0.)    ? 0. : 4. E+153;
  16.       v = n    = fint (m = x *    1.4426950406);
  17.       m -= n; if (*mx) ++(*mx); e =    chebgen    (m - 1., expcc,    ns);
  18.       v = (*ex += v); if (v    >= 0x400) return 4. E+153;
  19.       return (v <= 0) ? 0. : e; }
  20.  
  21. float loge (x)    float x;
  22.     { float     chebgen();
  23.       static float e, m, n,    logcc[]    =
  24.       {  9.3022922134 E-1,    -8.1841456644 E-2,   9.4766115856 E-3,
  25.         -1.2282837393 E-3,     1.6939529027 E-4,  -2.4301272033 E-5,
  26.          3.5827616763 E-6,    -5.3890849189 E-7,   8.2315218987 E-8,
  27.         -1.2726716189 E-8,     1.9871237790 E-9,  -3.1280003518 E-10,
  28.          4.9576937640 E-11 };
  29.       static int v,    *nx = &n, ns = sizeof (logcc) /    sizeof (float);
  30.       if (x    <= 0.) return -4. E+153;
  31.       n = x; v = *nx - 0x200; *nx =    0x200;
  32.       if (n    <= 0.8)    { --v; ++(*nx);    }
  33.       e = chebgen (n * 2.5 - 3.0, logcc, ns);
  34.       return e * (n    - 1.) +    (v * 0.69314718025); }
  35.  
  36. float tan (x)  float x;
  37.     { float    m, _stgen(), sqrt();
  38.       m = _stgen (x, 1); return m /    sqrt (1. - m * m); }
  39.  
  40. float cos (x) float x;
  41.     { float    _stgen(); return _stgen    (x + 90., 0); }
  42.  
  43. float sin (x)  float x;
  44.     { float    _stgen(); return _stgen    (x, 0);    }
  45.  
  46. float _stgen (x, t)  float x;
  47.     { float    fint(),    chebgen();
  48.       static float m, n, ccsin [] =
  49.       {  1.2762789622 E+0,    -1.4263078457 E-1,   4.5590080009 E-3,
  50.         -6.8293756748 E-5,     5.9248092843 E-7,  -3.3513958000 E-9,
  51.          1.3336392981 E-11 };
  52.       static int *nx = &n, ns = sizeof (ccsin) / sizeof (float);
  53.       n = x/360.; n    -= fint    (n + 0.5); if (*nx) *nx    += 2;
  54.       if ((*nx & 0x7FFF) > 0x200) n    = t ?
  55.           (n > 0. ?    -2. : 2.) + n :    (n > 0.    ? 2. : -2.) - n;
  56.       n *= (m = n);    if (*nx) ++(*nx); --n;
  57.       return m * chebgen (n, ccsin,    ns); }
  58.  
  59. float acos (x) float x;
  60.     { float    atan(),    sqrt();
  61.       return 90. - atan (x / sqrt (1. - x *    x)); }
  62.  
  63. float asin (x) float x;
  64.     { float    atan(),    sqrt();
  65.       return atan (x / sqrt    (1. - x    * x)); }
  66.  
  67. float atan (x)    float x;
  68.     { float    fabs(),    chebgen();
  69.       static float m, n, atancc[] =
  70.       {  8.8137358699 E-1,    -5.2946462263 E-2,   5.5679210277 E-3,
  71.         -6.9059750155 E-4,     9.2871486598 E-5,  -1.3107598052 E-5,
  72.          1.9105182955 E-6,    -2.8495930820 E-7,   4.3244389300 E-8,
  73.         -6.6516919840 E-9,     1.0342528811 E-9,  -1.6224319589 E-10,
  74.          2.5639831258 E-11 };
  75.       static int t,    *nx = &n, ns = sizeof (atancc) / sizeof    (float);
  76.       m = (t = (fabs (x) >=    1.)) ? -1./x : x;
  77.       n = m    * m; if    (*nx) ++(*nx);
  78.       n = m    * 57.295779513 * chebgen (n - 1., atancc, ns);
  79.       if (t) n = (x    > 0.) ?    n + 90.    : n - 90.; return n; }
  80.  
  81. float chebgen (zz, cc, n)  float zz, cc[];
  82.     { static float a0, a1, a2, *cp;    static int *ax = &a0;
  83.       cp = cc + n; a0 = a1 = 0.;
  84.       while    (cp > cc) { a2 = a1; a1    = a0; a0 = zz *    a1;
  85.                 if (*ax) ++(*ax); a0 += (*--cp - a2); }
  86.       return a0 - a2; }
  87.  
  88.  
  89. /* GENERAL FLOAT FUNCTIONS */
  90.  
  91. float fabs (b) float b;
  92. /* return absolute value */
  93.     { static unsigned *e; e    = &b;
  94.       *e &=    0x7FFF;    return b; }
  95.  
  96. int fsgn (b) float b;
  97. /* return signum: -1, 0    or 1 */
  98.     { static unsigned *e; e    = &b;
  99.       if (*e) return (*e & 0x8000) ? -1 : 1; return    0; }
  100.  
  101. float sqrt (x) float x;
  102. /* return square root */
  103.      {    if (x <= 0.) return 0.;    &x; inline (
  104.     0x89, 0xDF, 0x8B, 0x55,    0x04, 0x8B, 0x4D, 0x02,    0x8B, 0x1D, 0x31,
  105.     0xF6, 0x89, 0xD8, 0x81,    0xEB, 0xFF, 0x01, 0xD1,    0xFB, 0x81, 0xC3,
  106.     0x00, 0x02, 0x25, 0x01,    0x00, 0x74, 0x06, 0xD1,    0xE9, 0xD1, 0xDA,
  107.     0xD1, 0xDE, 0x56, 0x52,    0x51, 0x31, 0xC9, 0x31,    0xD2, 0xB8, 0x00,
  108.     0x40, 0x51, 0x51, 0x50,    0x89, 0xE5, 0xB8, 0x21,    0x00, 0x50, 0xD1,
  109.     0xD0, 0x35, 0x01, 0x00,    0x0B, 0x76, 0x04, 0x0B,    0x56, 0x02, 0x0B,
  110.     0x4E, 0x00, 0x29, 0x76,    0x0A, 0x19, 0x56, 0x08,    0x19, 0x4E, 0x06,
  111.     0x73, 0x02, 0xD1, 0xD8,    0x73, 0x14, 0x01, 0x76,    0x0A, 0x11, 0x56,
  112.     0x08, 0x11, 0x4E, 0x06,    0x33, 0x76, 0x04, 0x33,    0x56, 0x02, 0x33,
  113.     0x4E, 0x00, 0xEB, 0x09,    0x03, 0x76, 0x04, 0x13,    0x56, 0x02, 0x13,
  114.     0x4E, 0x00, 0xD1, 0x6E,    0x00, 0xD1, 0x5E, 0x02,    0xD1, 0x5E, 0x04,
  115.     0xD1, 0x66, 0x0A, 0xD1,    0x56, 0x08, 0xD1, 0x56,    0x06, 0x58, 0x48,
  116.     0x75, 0xAF, 0x83, 0xC4,    0x0C, 0x01, 0xF6, 0x73,    0x04, 0x42, 0x75,
  117.     0x01, 0x41, 0x89, 0x55,    0x04, 0x89, 0x4D, 0x02,    0x89, 0x1D ); }
  118.  
  119. float fint (b)    float b;
  120. /* return next int value below - eg -1.5 gives -2. */
  121.     { static float a; static unsigned aex;
  122.       a = b; aex = * (unsigned *) &a;
  123.       if ((aex & 0x7FFF) <=    0x200)
  124.       return (aex &    0x8000)    ? -1.0 : 0.0;
  125.       if ((aex & 0x7FFF) >=    0x220 )    return a;
  126.       /* shift down    mantissa by (0x220 - exp)
  127.          inc if negative and renormalise      */
  128.       inline ( 0xBD, &a,   0x31, 0xDB, 0xB9, 0x20, 0x02,
  129.            0x8B, 0x46, 0x00, 0x25, 0xFF, 0x7F, 0x29, 0xC1,
  130.            0xD1, 0x6E, 0x02, 0xD1, 0x5E, 0x04,
  131.            0x83, 0xD3, 0x00, 0xFF, 0x46, 0x00,
  132.            0x49, 0x75, 0xF1, 0x09, 0xDB, 0x74, 0x11,
  133.            0x8B, 0x5E, 0x00, 0x81, 0xE3, 0x00, 0x80, 0x74, 0x08,
  134.            0xFF, 0x46, 0x04, 0x75, 0x03, 0xFF, 0x46, 0x02,
  135.            0x8B, 0x46, 0x02, 0x0B, 0x46, 0x04,
  136.            0x75, 0x05, 0x89, 0x46, 0x00, 0xEB, 0x12,
  137.            0xF7, 0x46, 0x02, 0x00, 0x80, 0x75, 0x0B, 0xD1, 0x66,
  138.            0x04, 0xD1, 0x56, 0x02, 0xFF, 0x4E, 0x00, 0xEB, 0xEE    );
  139.       return a; }
  140.  
  141.  
  142. /* Convert string representation to float */
  143. float atof (s) char *s;
  144.      {    static int e, m; float f, g, fe10();
  145.     m = e =    0; f = 0.; g = 1.;
  146.     while (*s == ' ' || *s == '\t')    ++s;
  147.     if (*s == '-' || *s == '+') m =    (*s++ == '-');
  148.     if (*s && *s !=    '.')
  149.        while (*s >=    '0' && *s <= '9')
  150.        f = f * 10. + (*s++ - '0');
  151.     if (*s == '.') { ++s;
  152.        while (*s >=    '0' && *s <= '9')
  153.        f +=    (g *= .1) * (*s++ - '0');  }
  154.     if (m) f = -f; m = 0;
  155.     while (*s == ' ' || *s == '\t')    ++s;
  156.     if (*s == 'E' || *s == 'e')
  157.     {  ++s;    while (*s == ' ' || *s == '\t')    ++s;
  158.        if (*s == '-' || *s == '+') m = (*s++ == '-');
  159.        while (*s >=    '0' && *s <= '9')
  160.          e = e * 10    + (*s++    - '0');
  161.        if (m) e = -e; f = fe10 (f, e); }
  162.     return f; }
  163.  
  164.  
  165. /* PRINTF/SPRINTF/FPRINTF ROUTINES - for FP print capture */
  166. printf (fs)int *fs; {int _psc(), _fprt(); return _prout(stdout,    fs, _psc); }
  167. sprintf(fs)int *fs; {int _pss(), _fprt(); return _prout(fs,   --fs, _pss); }
  168. fprintf(fs)int *fs; {int putc(), _fprt(); return _prout(*fs,  --fs, putc); }
  169.  
  170. /* Floating point routine for _prout() */
  171. _fprt (f, b, p,    t) float f; char *b;
  172.      {    if (t) p *= t; else p =    (f < 1000000.) ? p : -p;
  173.     if ((p = ftoa (f, b, p)) < 0 ||    t) return;
  174.     b += p;    while (*--b == '0') *b = ' '; }
  175.  
  176.  
  177. /* FP to Ascii representation.    Exponential format used
  178.    if precision    is negative, or    abs (f)    exceeds    1. E 9.
  179.    Returns +/- string length, guaranteed less than 30 chars. */
  180. ftoa (f, s, p) float f;    char *s;
  181.      {    static int dd[3], ee, ef, sg, zz, qq;
  182.     static char *t;    float fe10();
  183.     struct { unsigned _exp,    _mhi, _mlo; };
  184.  
  185.     if ((qq    = p) < 0) { ef = 1; p =    -p; } else ef =    0;
  186.     if (p >    14) p =    14; sg = f._exp    & 0x8000;
  187.     if ((f._exp &= 0x7FFF) > 0x21E)    ef = 1;
  188.  
  189.     if (f._exp)
  190.     {  if (!ef) f += fe10 (.5, -p);
  191.        ee =    f._exp - 0x202;
  192.        ee =    (3 * ee    + ee / 100) / 10 + (ee > 0 ? 1 : 0);
  193.        f = fe10 (f,    -ee);
  194.        if (f._exp >= 0x201)    { f = f    * .1; ++ee; }
  195.        if (ef)
  196.        {  if (qq >=    0) p +=    2; f +=    fe10 (.05, -p);
  197.           if (f._exp >= 0x201) { f = f * .1; ++ee; } }
  198.        dd[0] = f._mhi; dd[1] = f._mlo; dd[2] = 0;
  199.        if ((zz = 0x200 - f._exp) > 0)
  200.           while (zz--) /* shift array right    */
  201.           inline ( 0xBB, dd,   0xD1, 0x2F, 0x43, 0x43,
  202.                0xD1, 0x1F, 0x43, 0x43, 0xD1, 0x1F ); }
  203.     else { ee = 1; dd[0] = dd[1] = dd[2] = 0; }
  204.  
  205.     t = s +    1; zz =    0;
  206.     if (!ef) { if (ee <= 0)    *t++ = '0';
  207.            else    do *t++    = '0' +    _fpc (dd, &zz);    while (--ee); }
  208.     else *t++ = '0'    + _fpc (dd, &zz);      *t++ = '.';
  209.     if (!ef    && ee) while (ee++ && p-- > 0) *t++ = '0';
  210.     while (p-- > 0)    *t++ = '0' + _fpc (dd, &zz);
  211.     if (ef)     { *t++    = ' '; *t++ = 'E';
  212.            if (--ee >= 0) *t++ = '+';
  213.            else    { *t++ = '-'; ee = -ee;    }
  214.            if (ee >= 100) { *t++ = '0' + ee/100; ee %= 100; }
  215.            *t++    = '0' +    ee/10; *t++ = '0' + ee % 10; }
  216.     *s = (zz && sg)    ? '-' :    ' '; *t    = 0;
  217.     return ef ? s -    t : t -    s; }
  218.  
  219. /* Unload characters from array    - or with c */
  220. _fpc (a, c) int    *a, *c;
  221.      {    inline ( 0x89, 0xE5, 0x8B, 0x76, 0x04, 0x83, 0xC6, 0x06,
  222.          0x31, 0xDB, 0xB9, 0x03, 0x00, 0x4E, 0x4E, 0xB8,
  223.          0x0A, 0x00, 0xF7, 0x24, 0x01, 0xD8, 0x89, 0x04,
  224.          0xBB, 0x00, 0x00, 0x11, 0xD3, 0xE2, 0xEE,
  225.          0x8B, 0x76, 0x02, 0x09, 0x1C ); }
  226.  
  227. float fe10 (b, p) float    b;
  228. /* return b * 1. E+p */
  229.     { static float m;
  230.       if (p    > 0) m = 10.; else { p = -p; m = .1; }
  231.       while    (p) { if (p & 0x1) b *=    m;
  232.               if (p >>=    1) m *=    m; } return b; }
  233.  
  234. /* End of float    library    */
  235. #list+
  236.